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For the kinetic energy of Id model finite systems the leading corrections to local approximations 
as a functional of the potential are derived using semiclassical methods. The corrections are simple, 
non-local functionals of the potential. Turning points produce quantum oscillations leading to energy 
corrections, which are completely different from the gradient corrections that occur in bulk systems 
with slowly- varying densities. Approximations that include quantum corrections are typically much 
more accurate than their local analogs. The consequences for density functional theory are discussed. 



I. INTRODUCTION 

Modern density functional theory (DFT) has become 
a popular electronic structure method, because of its bal- 
ance between computational efficiency and accuracy. The 
original density functional theory was that of Thomas 1 
and Fermi 2 (TF), in which all parts of the electronic 
Hamiltonian are approximated by explicit density func- 
tionals, and the energy minimized over possible densities. 
Typical errors in TF theory are of order 10% of the total 
energy, but molecules don't bind. 3 In the 1950's, Slater 4 
introduced the idea of orbitals in DFT, i.e., solving a set 
of single-particle equations to construct the energy of the 
interacting system, which is typically much more accu- 
rate. This was shown to be a formally exact approach in 
the celebrated works of Hohcnberg and Kohn 5 and Kohn 
and Sham. 6 The latter also introduced the local den- 
sity approximation (LDA) for the only unknown needed 
to solve the Kohn-Sham (KS) equations, the exchange- 
correlation (XC) energy as a functional of the density. 
Since then, the field has gradually evolved with improve- 
ments in computational power, algorithms, and approxi- 
mate functionals to the workhorse it is today. 7 

Unfortunately, the existence theorems give no hint of 
how to construct approximate functionals. Presently, 
there is a dazzling number of such approximate func- 
tionals suggested in the literature, and implemented in 
standard codes, both in physics and chemistry. 8 Many 
of these are physically motivated, and work well for the 
systems and properties for which they were designed, but 
usually fail elsewhere. There appears to be no systematic 
approach to the construction of these functionals, beyond 
artful constraint satisfaction. 9 

In the present paper, we return to the origins of DFT 
and ask, what are the leading corrections to the origi- 
nal approximation of Kohn and Sham, the LDA? This is 
a very difficult question that we can only aspire to an- 
swer for the XC energy for any electronic system. In the 
present article, we answer the question for an extremely 
simple case, but one that contains many features relevant 
to the problems of electronic structure. 

The original works on density functional approxima- 
tions emphasize the gradient expansion, 5,6 which is a par- 
ticular approach to improve upon a local density approx- 
imation. Imagine an infinitely extended slowly-varying 



gas. The corrections to the local approximation are given 
accurately by leading corrections in the density gradi- 
ent. But real systems do not look like slowly- varying 
gases. All finite systems have evanescent regions, as do 
many bulk solids. The regions can be separated via clas- 
sical turning-point surfaces, evaluated at the Fermi en- 
ergy of the system. 10 Typically, such regions are atomic 
sized. Most importantly, the gradient expansion fails 
completely both near and outside these surfaces. Gener- 
alized gradient approximations (GGAs) and other meth- 
ods have been developed to overcome these difficulties. 
These include only a finite order of gradients, but employ 
a form which contains many powers those gradients. 

To study the effects of confinement to finite re- 
gions on density functional approximations, we use non- 
interacting particles, and study only their kinetic energy, 
which was locally approximated in the original TF the- 
ory. We study only one dimension, where semiclassi- 
cal approximations are simple, and the WKB 11-13 form 
applies in the absence of classical turning points where 
the potential v(x) has finite slope. We avoid such turn- 
ing points by using box boundary conditions and study- 
ing only systems whose chemical potential is everywhere 
above v(x). 

The answer is surprising: For most systems, the lead- 
ing corrections (in a sense that will become clear within) 
are not the simple gradient corrections commonly dis- 
cussed, and used as starting points to construct GGAs. 
Instead, both the density and kinetic energy density can 
be very accurately approximated as functionals of the 
potential. 

The limit we discuss is also an important result in it- 
self. We carefully show precisely how TF becomes exact 
in a semiclassical limit. Essentially, we take H — > 0, keep- 
ing the chemical potential /i roughly fixed, and renormal- 
izing the density so as to retain the original particle num- 
ber. If, further, one performs a moving-average over the 
space coordinate, with a range chosen to be small com- 
pared to the spatial variation of the potential, but large 
compared to quantum oscillations as h — > 0, the density 
uniformly converges to that of TF theory. We call this 
the continuum limit of a finite system. The separation 
between quantum eigenvalues becomes infinitesimal, and 
all sums become integrals. The integrands within con- 
tain purely classical quantities in terms of the potential, 
v(x). A similar simplification occurs for the kinetic en- 
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ergy density, given in terms of v(x), and when v(x) is 
eliminated from the two expressions, what remains is the 
LDA to the kinetic energy. 

Having carefully defined this limit, we can then discuss 
the approach to that limit, and the leading corrections 
to the local approximation. We find that the dominant 
corrections (in 1/h) are not gradient corrections due to 
the variations in v(x) in the interior, but rather are quan- 
tum oscillations due to the hard walls at the boundaries. 
These quantum oscillations are generic features of any 
quantum system, and their nature is determined by the 
classical turning points. 14 They give rise to the phase cor- 
rections to the classical density of states in the Gutzwiller 
trace formula, 15 and will be present for any finite quan- 
tum system. The only case in which they vanish is that 
of periodic boundary conditions with the chemical po- 
tential above the maximum of the potential. Only in 
such systems does the gradient expansion produce the 
correct asymptotic expansion in powers of H, equivalent 
to gradients of the potential. For any finite system, the 
series eventually diverges, but truncation at a lower or- 
der can yield highly accurate results, if the gradients are 
sufficiently small. 
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FIG. 1: Exact, TF, our semiclassical (sc), and WKB density 
for a single particle in v(x) — — 12 sin 2 (7ra;) with hard walls at 
x = and x = 1 (h = m = 1). 

To give some idea of the power of the methods we 
develop, in Fig. 1 we show the density of one particle 
in a simple well (v(x) = — Dsin 2 (nx/L) in a box from 
to L, with D = 12.) The exact density is found by 
numerically solving the Schrddinger equation. The TF 
density is found by minimizing the Id TF kinetic energy 
functional and choosing the chemical potential to yield 
one particle. That density smoothly follows the potential. 
Due to the hard walls there are no classical turning points 
where v(x) has finite slope. Hence, a WKB treatment 
can be applied here, yielding a WKB eigenvalue that is 
positive and reasonably accurate. But the result of our 
present analysis is a simple formula for the density, which 
uses WKB wavef unctions as input, is much more accurate 



still. Perhaps more importantly, we have been unable 
to create situations where our approximation completely 
fails. 

In the electronic structure problem, the local approxi- 
mation to the XC energy is analogous to the TF approx- 
imation in Fig. 1. While there are many excellent ap- 
proximations that improve over the local approximation 
by typically an order of magnitude, they are usually tai- 
lored to specific systems or properties, and contain either 
empirical parameters 16 or at least careful selection of ex- 
act conditions to impose on an approximation to fix the 
parameters non-empirically. 17 Our formulas are derived 
via a semiclassical analysis that yields unique approxima- 
tions, are robust, and typically two orders of magnitude 
better than the local approximation. 

The semiclassical result for the density was given in 
a short report 18 , but the kinetic energy density formula 
derived there failed close to the boundaries. Here we ex- 
plain that failure, in terms of boundary-layer theory 19 ' 20 , 
but applied to sums over eigenstates rather than individ- 
ual solutions to a differential equation. By identifying 
our limit, we can then cleanly separate the two differ- 
ent length scales in the problem. Earlier approaches 14 ' 21 
yield only the asymptotic result in the interior of the 
system (ironically labelled the outer region in boundary 
layer theory), but that there exists a region close to each 
wall (appropriately called a boundary layer) where that 
solution fails. However, within the boundary layer, a 
different asymptotic expansion applies and, by matching 
these two solutions, we construct a uniform approxima- 
tion that is asymptotic to a given order everywhere in 
the system. Many different aspects of these issues have 
been addressed over the decades since the original work of 
Kohn and Sham. 14 For example, Balian and Bloch, 22 in 
the context of nuclear physics, identified the need for spa- 
tial averaging to approach the limit. In the early 1970's, 
Yuan and Light 23,24 and coworkers 25 developed the the- 
ory in terms of path integrals and density matrices. Re- 
cently, Ullmo, Barangcr, and coworkers 26,27 studied the 
nature of quantum oscillations in application to quantum 
dots. 

What is the significance of our results for the real world 
of 3d, interacting electrons? Our results, for a very dif- 
ferent case, reveal the nature of the corrections to local 
approximations. These will differ in detail depending on 
the dimensionality of the system, or the specific func- 
tional being approximated. However, qualitative features 
(such as the local approximation becoming exact in the 
classical continuum limit, gradient expansions being in- 
valid near turning points, etc.) are general. Thus our 
analysis can (and already has 28 ) provided guidance for 
the construction of new XC density functionals. On the 
other hand, there is also a considerable amount of work 
done in the field of orbital-free DFT, 29 ' 30 but effort is fo- 
cused on finding an accurate approximation to the non- 
interacting kinetic energy as an explicit functional of the 
density. The present work derives potential functional 
approximations, an entirely different matter, and so has 
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no overlap with existing work in that field. Our work 
suggests that the potential is a better variable than the 
density, and we show how corrections to local approxi- 
mations of the density and kinetic energy density can be 
derived as potential functionals for simple model systems, 
but the methods and results shown here do not readily 
generalize to three dimensions. 

This paper is divided as follows: In Sec. II, we intro- 
duce our notation and define the continuum limit and 
show that local approximations become exact in this 
limit. Next, in Sec. Ill we derive the leading corrections 
for both the density and kinetic energy density as func- 
tionals of the potential by explicit summation of WKB 
orbitals. Then, in Sec. IV we "fix" the difficulties at the 
walls to produce a uniform approximation everywhere, 
and then study its properties comparing to the exact 
result both in the classical continuum and the large-A 
limits. Finally, in Sec. V we discuss the implications and 
relevance for real electronic structure calculations. In 
the appendix we give a detailed derivation of the inte- 
rior solution of the density and the KED using the WKB 
Green's function in the complex plane, just as has been 
done before. 



II. CLASSICAL CONTINUUM 

In this section, we introduce our notation and briefly 
review the salient points known from the literature. 
As discussed in the introduction, we restrict ourselves 
to non-interacting particles in one dimension. We use 
atomic units throughout (e 2 = h = m e = 1), so that 
lengths are expressed in Bohr radii, and energies in 
hartree. 



A. Background and notation 

We write the Hamiltonian as 



1 d 2 



(1) 



where a hat denotes an operator with t being the kinetic 
energy operator, and v(x) the potential. We denote the 
solutions to the Schrddinger equation as 



h(j>j(x) = ej (j)j{x), i = l,2, 



(2) 



The solutions to the Schrodinger equation can be ex- 
panded in powers of fi, 31 and retaining just the first two, 
we find the WKB solutions for a given energy e are 



/WKB/ \ 

9 \x) 



i 8(x) 



(3) 



and its complex conjugate, where the dependence on e is 
via the definitions of the wavevector 



k(x) = ^2(e-v{x)) 



(4) 



and classical phase 



6{x) 



dx' k(x') 



(5) 



where the constant is arbitrary. These solutions are exact 
when the potential is constant, and highly accurate when 
the potential is slowly varying on a scale determined by 
the energy. However, the particular combination that 
forms an eigenstate depends on the boundary conditions. 
The density is then 



n(*)«-/ &|0 WKB ( £ ,z)| 2 = ^ 



1 

n 



(6) 



where k^x) is the wavevector evaluated at /z, the chem- 
ical potential for the system, found via normalization. 
Note that k^{x) is a function of /i — v(x) alone, so we 
define the local chemical potential: 



fx(x) = IX — v(x). 



(7) 



Then, because WKB is exact for an infinitely extended 
system with constant potential (free particles), we find 
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n™*(n) = -(2 M )^, t™«(») = ^{2^f 2 . (8) 

7T 07T 

The corresponding integrals are local potential approxi- 
mations (LPAs) to the particle number and kinetic en- 
ergy: 



N loc [n{x)} = J dxn UDii {n{x)), 
T loc [fx(x)] = [ dxt unil {^{x)). 



(9) 



Inverting the relation n u [fx(x)] and inserting into 
i unlf |/t(x)] yields the local density approximation to T: 



T lo >]= dxt unit (n{x)), t unii (n(x)) = 



tv n d (x) 
6 



(10) 



This is the one-dimensional analog of the TF kinetic en- 
ergy density functional (up to simple factors of 2 for dou- 
ble occupation). 32 

One can also work backwards from the LDA to the 
WKB results. The LDA for the kinetic energy allows us 
to find an approximate density for a given potential, by 
minimizing the total energy subject to the constraint of 
a given particle number. This produces 



-n 2 (x) + v(x) = fi . 



(11) 



the TF equation for this problem, which is identical 
to Eq. (6). The solution is the TF density, n TF (x) = 
n umf (fj,(x)). The total particle number A is a continuous 
monotonic function of the parameter \i that is invcrtiblc 
for \i > v m i n , the minimum of the potential. 



4 



The leading corrections to the WKB wavefunctions, 
i.e., the next two powers in ft, are well-known 19 and 
produce constant corrections to both the phase and the 
wavevector. Samaj and Percus 32 showed very elegantly 
how the series can be generated to any desired order. 
Continuing with the higher-order corrections to WKB 
leads to corrections that depend on derivatives of the po- 
tential, where v'(x) = dv/dx. The potential gradient 
expansion for the density is 



>(x)]( a 



\n{x)) 1 - 



v"(x) v' 2 (x) 



+ 



and for the kinetic energy density 



(12) 



^ (a5 ) ] („ = ^) ) (i-^ + ^M + . 

(13) 

which, when inverted leads to the density gradient ex- 
pansion for T: 

T[n] w y J dxn 3 (x)--^ J dx 



by 7ft. As 7 becomes very small, the discrete levels of the 
system merge, and the envelope of their density of states 
approaches a well-defined limit. We call this limit the 
classical continuum. While it has been long understood 
that local density approximations become exact in this 
limit, 35 relatively little attention has been paid to how 
exactly this limit is reached in a quantum system. 

Consider a Id box of length L with given potential v(x) 
and particle number N, i.e., the lowest TV eigenstates are 
occupied. Then increase the particle number to N', but 
\ choose 

' 7 = N/N 1 < 1 (15) 

i.e., ft is reduced in proportion to the increase in particle 
number. Of course, there will now be N' particles in our 
well, so define 



\x) 



2n(x) 



'(*) 



(14) 

A gradient expansion approximation is the finite trunca- 
tion of that series. Because the semiclassical expansion is 
asymptotic, this is an asymptotic expansion for periodic 
systems where fj, is above the maximum of the poten- 
tial. It can be made arbitrarily accurate by application 
to sufficiently smooth densities, but for any given den- 
sity, addition of sufficient terms will eventually lead to 
divergence. For example, a potential that contains steps 
will produce divergences beyond the lowest order. Also, 
Coulomb potentials are known to vary too rapidly for 
such expansions to apply. 33 We note that, in Id, because 
of the negative coefficient in the gradient correction, min- 
imizing the total energy is unbounded and nonsensical in 
the presence of this correction. 34 We also note that ft 
never appears in the functional dependence on the den- 
sity in Eq. (14). 



- ( \ N r ^ 



(16) 



as a renormalized density, whose particle number 
matches the original value at 7 = 1. This process is 
illustrated in Fig 2, where we plot renormalized densi- 
ties for several particle numbers N', and the TF result 
in the same potential as used for Fig. 1. One can see 
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B. The classical continuum limit 

We define a continuum as any region of energy in which 
the eigenvalues of the Hamiltonian are not discrete. The 
first, simplest example is that of a particle in a well, with 
a potential set that vanishes as |ar| — > 00. For e > 0, there 
is the free-particle continuum, with scattering states of 
the system that cannot be box-normalized. Another con- 
tinuum arises in solid-state physics, when we apply pe- 
riodic boundary conditions to our potentials, in order to 
simulate bulk matter. Then, for single-particle states, the 
energy levels form distinct bands, usually labelled by a 
wavevector. Within each band, the energy is continuous. 
We call this the bulk or thermodynamic continuum. 

But any system also has a classical continuum, which 
can be found by letting 7 — > 0, where we have replaced ft 



FIG. 2: TF and renormalized exact densities for N' = 
1, 4, and 16 particles in v (x) = -12 sin 2 (the), < x < 1, 
showing approach to continuum limit. 

how our procedure reproduces the TF density, almost. 
As N' grows, the oscillations in the interior of the box 
become smaller (with an amplitude of 0(1/N')), while at 
the edges (within 0(L/N') of the wall), the exact density 
always drops to satisfy the boundary condition, while the 
TF density does not. So we also define a moving average 
of a function of a; as: 

m+Ai/2 

(n(x)) A x = / dx'n{x')/Ax. (17) 

Jx-Ax/2 

The length scale of the moving average is chosen to be 
much larger than that of the quantum oscillations of the 
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exact density and of the boundary region at the wall, but 
still vanishes as 7 — > 0. Then, finally, 

lim(n 7(JV '^)(M,a;)>V7Z, = njf(x). (18) 

Thus we see that the TF density in a given problem 
is the limit as 7 — > 0, but the convergence is highly non- 
uniform. At the walls, the true density is always zero, 
but the TF density is finite. There are likely many other 
averaging procedures, such as taking the limit as a finite 
temperature vanishes, which can be used to define the 
limit, but the current one is suffificicnt for our present 
purpose. Similarly, we define t 7 (x) = 7 3 i 7 (x). 



expansion in h. We wish to develop approximations that 
are smooth in ft, and yield the exact approach to the clas- 
sical continuum limit, ignoring the discrete nature of the 
eigenstates, i.e., we wish to build in the smooth envelope 
of functions such as the density. At the end, we requan- 
tize our results, and find more accurate summations than 
WKB, even for N = 1. 

Begin with the smooth (non-oscillating) contribution. 
We use the Euler-Maclaurin formula in the following 
form: 



N 



djfj - -^(/ F 



£) + O(f') , (22) 



III. LEADING CORRECTIONS TO LOCAL 
APPROXIMATIONS 

In this section, we derive the leading corrections in 
7 scmiclassically, using only elementary techniques, for 
the sake of transparency. The first such derivation was 
by Kohn and Sham, 14 using a very elegant analysis of 
the properties of the Green's function in the complex 
plane. We include an appendix in which we also derive 
our results via this method. In this section, we simply 
derive formulas for large N and extract the 7-dependcncc 
from such formulas. 



A. Density 

As in Sec. II, the density of N non-interacting fermions 
is approximated by the sum of the squares of the WKB 
orbitals, normalized and satisfying the boundary condi- 
tions. Because we are deriving the leading corrections, 
we carefully normalize here: 



sin 6j (x) , 



(19) 



where the normalization constant is found by ignoring 
the oscillating term. Define 



Tj(x) 



dx 
(x) 



T 3 =t {L), 



(20) 



the classical time for a particle in level j to travel from 
to L. Thus, in WKB theory, the density is approximated 
by 



N 



WKB 



1 — cos29j(x) 



E kj(x)Tj 



(21) 



Performing the sum exactly yields nothing other than the 
standard WKB approximation to the density as the sum 
of WKB densities 36 . However, such an approximation is 
inconsistent, since it retains the discrete nature of the 
eigenvalues, and will not yield a sum with a well-defined 



where prime denotes a derivative with respect to j and 
a subscript F denotes evaluation at the upper limit of 
the integral, j F = N + 1/2, while a subscript m denotes 
evaluation at the lower limit, j = 1/2. This is an ex- 
pansion for sums in the same parameter as for the WKB 
cigenfunctions, i.e., gradients of the potential. We re- 
tain only the first two terms, consistent with our WKB 
approximation for the orbitals. 

To expand the sum of the smooth terms in such pow- 
ers, we need to relate the level index j with the energy 
in a continuous fashion. Write the WKB quantization 
condition as 



® j =6{e jt L)=jn, 



5 = 1,2, 



(23) 



which defines ij, the WKB eigenvalue implicitly. Then 
differentiation yields: 



e j = T i l 'i = 71 ■ 



(24) 



where fx sc = e F = £n+i/2- This allows us to apply 
Eq. (22) to the smooth contribution from Eq. (21). De- 
fine 



Kj(x) 



1 



kj ( x ) Tj 



(25) 



which has units of inverse length and whose j-dependence 
is typically weak, vanishing entirely for a flat box. Then 



N 

E K j(- T ) 



djKj = 



dk 

km * 



(26) 



where we have neglected terms that contain derivatives 
of Kj at the end points. Thus 



Kj(x) w (k F (x) — k m (x)) /it , 

3=1 



(27) 



where the quantities in the integrand of Eq. (26) depend 
on j in a continuous manner, k F (x) = ^/2/i sc (x) and 
fi sc satisfies the quantization condition in Eq. (23) with 
j = N + 1/2, while k m (x) is the same but with N = 0. 
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The oscillating term in Eq. (21) is more delicate. For 
each x, we expand 6j (x) about its value at the Fermi level 
linearly 



9j(x) = e F (x) - (j F - j) a F (x) + ..., (28) 



where 



otj{x) = O'Ax) = 7T 



(29) 



If we truncate at this level, and use the geometric sum 
defined by 



/v 



1-z 



Z J = z ■ 



1 - z 



(30) 



and using z = exp[i 2a F (x)}, we find 



N 

E 



cos29j(x) k m (x) K F (x)sm20 F (x) 



2 sina F (a;) 



(31) 



The first term here exactly cancels the second term of 
the smooth contribution in Eq. (27). It is found from 
performing the geometric sum, using Eq. (28) to undo 
the linear approximation at the lower end of the sum in 
Eq. (31). Hence, the semiclassical density is 



n sc (x) = n s (x) + n c 
where s denotes the smooth term, 



n s (x) 



k F (x) 



(32) 



(33) 



and osc the oscillating contribution, defined to have zero 
moving average as 7 — > 0, 



=(X) 



sin2# F (:r) 



2T F k F (x) sin a F (x) 



(34) 



Note how completely different these corrections are 
from those of Eq. (12). This is a central result of this 
work. 

In general, the smooth term does not match that of 
TF theory, because it is evaluated at N + 1/2, not N. 
This is not an artifact of Eq. (22), but reflects the 1/2 
electron loss of density in the quantum correction. We 
write Eq. (22) in a form that avoids terms with f m and 
/ F , but the 1/2 term is independent of any particular 
choice. 

Note also that our semiclassical density is not normal- 
ized, in general. For a flat box, the quantization does 
imply correct normalization, but not more generally. It 
is straightforward to find slightly modified definitions of 
0(x), etc., that both normalize the density and satisfy 
the boundary conditions (i.e., 9 = jn), but we choose to 
retain this error as a measure of the error in our semi- 
classical approximations. We discuss this fact and assess 
the error in Sec. IV A. 



Finally, we can rewrite the result as: 

k F (x) 



n sc (x) 



[l-r)(x)f(a F )w(e F )] 



where f(a) = 1/sina, w = sin(2#), and 
7r huj F 



r)(x) 



2k 2 (x)T F 8(/i sc - v{x)) 



(35) 



(36) 



is the small parameter, once x is not too close to the 
wall, and uj f = 2-k /T F is the classical frequency of colli- 
sions with the walls at the Fermi energy. To show the 
7-dependencc explicitly, replace N by N/"f and write 



n sc , 7 (a;) = jn sc , K (x) 



(37) 



7 



sin26' F!7 (a;) 



7r 2T F:7 k Fn (x) sin a F ^(x) 
where F, 7 denotes evaluation at N/j + 1/2. 

B. Kinetic energy density 



A similar analysis can be applied to the kinetic energy 
density, but must be done more carefully: 

*(*) = E ^rV-(z)! 2 « E - ^(z)) , 



(38) 

where £j(x) = kj(x) Kj(x). First we evaluate the sum 
over the smooth contribution using the same logic as for 
the smooth piece of the density. Applying Eq. (22) we 
obtain 



N 

E 

i=i 



6tt 48 



0(0 



(39) 



We know that the contributions from the lower end will 
be cancelled by analogous contributions in the oscillating 
piece. To evaluate that, we define: 



N 



h (p) {z) = EO'f - Jf z j = z hW{z) - j F h<?-V{z) , 
i=i 

(40) 

where h^'(z) = dh^/dz. 

Each term has many terms, but only those containing a 
z N will contribute to our answer, because when we insert 
Eq. (28) into Eq. (38), the prefactor contains z~ N . Then 



2i{e F -j F a F ) 



tosc (*e) 

Evaluating term by term yields 
iosc(z) = 77 Uf / 



^(°>-a (1) + ^ (2) 
(41) 
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df dw 1 9 2 / 



da 89 



da 2 



(42) 
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The derivatives of £ w.r.t. j can be written as 
and 



2.„// 
3 



(43) 



(44) 



boundary conditions will always be satisfied, no matter 
what v(x) is. Outside that region, i.e., in the interior of 
the box [xp < x < L — xp), the nonlocal k F (x) is used. 
We illustrate our approximations and the exact KED for 



We now drop all derivatives of Kj, because they vanish 
in the fiat box limit, yielding: 



t BC (x) = £ 



i 



3 i 
4^- 



d£&w _ 
da 86 



(45) 



Again, the explicit 7-dependence of this formula is found 
by replacing N by N + 1/2. As we shall show in later 
sections, this result is less well-behaved than that for the 
density. For example, when the potential is non-uniform, 
the semiclassical kinetic energy density of Eq. (45) incor- 
rectly fails to vanish at the edges. 

To overcome this failure, we define the edge as being 
those values of x up to some fraction (3 of a period of the 
classical phase: 



>(xp) = /3w, 



(46) 



such that the edge region is x < xp and x > L — xp. 
We choose (3 = 1/4 (and the interior is all the rest). 
This mimics the approach used in boundary-layer the- 
ory for differential equations, which can be applied to 
the fi-expansion of the individual levels. One constructs 
approximations that are correct to a given order in 
the asymptotic expansion in each region separately, and 
hopes to find a middle region where they match, yield- 
ing a solution with uniform convergence properties, 19 i.e., 
with the correct asymptotic expansion for any x. The 
only difference here is that we are applying these ideas 
to the sum of levels, not the individual levels themselves. 
Our final semiclassical approximation for the KED is 



tsc(x) 



^SC (^) 



if 

else , 



xp < x < L — xp . 



(47) 



where 



C if (*) = 



(fc™ lf ) 3 (fcJJ nlt ) 2 sin(2fc™ It a;) 

6 7r 4Lsin(7ra;/L) 
7T fc™ if cos(7ra;/i) cos(2 k^ nif x) 

4L 2 sin 2 (7ra/L) 
7T 2 sin(2 k™ if x) (\ 1 
8L 3 sin(7ra;/L) \2 ~ sm 2 {7Tx/L) 



.(48) 



Hence, inside the edge region we approximate the KED 
by i" c n (x) meaning that we evaluate Eq. (45) with the 
local Fermi wavevector for a uniform potential, i.e., re- 
placing k F everywhere by k F nil = y/2[i sc and defining all 
other quantities based on that. In particular the classi- 
cal phase and transit time become linear in s, as k F 
for the uniform system is independent of x. Hence, the 




^8 



FIG 

v(x) = — 12 sin 2 (7rx), where < x < 1 
value is eo = —4.27 and ^t sc = 5.52. 



3: Exact and approximate ground-state KEDs for 

The lowest eigen- 



a single- well potential v(x) — — 12 sin 2 (772;) within box 
boundaries in Fig. 3. Note that our present approxima- 
tion of Eq. (47) is substantially more accurate for both 
the edge region and the interior than the previously de- 
rived KED of Ref. 18. In the next section we discuss this 
fact more quantitatively. 



IV. PROPERTIES 

We next test our approximations, to demonstrate both 
their accuracy and that they have the properties claimed 
for them. We begin with several integrated quantities, 
mostly energies. 

A. Energies and normalization 

The LPA yields densities that are local in the potential, 
with the exception of the value of the chemical potential, 
which must be determined globally. The quantum cor- 
rections depend on several other terms, such as 9 F (x) and 
T F , which are still simple functionals of the potential, but 
distinctly non-local, depending on integrals over v(x). To 
test the integrated properties of the density, we calculate 
moments over that density. The obvious one is the third 
moment, as that is simply related to the local density 
approximation to T, evaluated on that density. 

We choose a standard potential, v(x) = —10 sin 2 nx in 
a box of length 1, and insert one particle. Both exact 
and approximate results are given in Table I. 

First note that the TF result, T TF , is about 50% too 
small, compared to the exact answer, T. This is the result 
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TABLE I: Exact and approximate quantities for one particle 
in a single-well potential v(x) = —10 sin 2 (ttx), < x < 1. T 
is the exact kinetic energy and n the exact density. 



Energy 


levels 


ei 


£2 


M 


A*sc 






-2.71 


14.6 


0.637 


6.38 




Kinetic 


energy 


T 


T loc [n] 






exact 


sc 


TF 


exact 


sc 






5.07 


5.02 


2.31 


4.93 


5.07 



of minimizing the energy using LDA as in the first term 
of Eq. (14). 

We measure the quality of the TF density and our 
semiclassical density by evaluating the LDA kinetic en- 
ergy on those densities, i.e., T loc \n TF ] = T TF and 




rplt 



J, where the point of reference is the LDA ki 



netic energy evaluated on the exact density, T loc [nl. But FIG - 4: LDA kinetic ener g v multiplied by the scale factor 
the TF result remains about 50% too small compared X for diff f e , nt 7 evaluated on n, (*), our semiclassical den 



to r u 

sity, T 



locr 



However, the LDA on our semiclassical den- 
i sc ], yields an energy only 3% too large, i.e., 
reducing the error by about a factor of 20. 

To test our semiclassical kinetic energy, T sc , we com- 
pare with the exact value, T, and find an error of only 
0.9% too small, i.e., 50 times better than T TF . Thus 
the semiclassical results are more than an order of mag- 
nitude better than bare DFT results, because they in- 
clude quantum oscillations. In fact, the LDA kinetic en- 
ergy evaluated on the exact density, T loc [n] , yields only a 
2.7% underestimate, showing that local approximations 
do much better on accurate densities, but still not as well 
as our direct approximation, T sc . 

These systems do not appear to be particularly semi- 
classical: the potential is not flat nor is the particle num- 
ber or index high. We can analyze the source of this 
accuracy by expanding integrated quantities in powers 
of 7 about 0: 

T( 7 ) = T<°> + 7 T« + 7 2 T< 3 > + ... (49) 

For the kinetic energy, from the previous discussion: 

T (0) = T TF (50) 

while our derivation should yield: 

T« = T« . (51) 

These results should hold for both the local approxima- 
tion applied to the exact density (and so test our semi- 
classical density) and the exact kinetic energy (and so 
test our semiclassical kinetic energy density). In Fig. 4, 
we study the 7-dependence of T loc [n] applied to various 
densities for a generic well. Clearly TF gives the 7 = 
value, while the semiclassical density includes the correct 
linear contribution, and is quite accurate for higher-order 
contributions. We also note that inclusion of the linear 
term greatly improves over the TF result, but that the 
LDA kinetic energy evaluated on our semiclassical den- 
sity is even more accurate still. 



sity rise, 7(2;), and the exact density for a single- well potential 
v(x) = —10 sin 2 (ttx). 



Because our expansion is in powers of h, we expect that 
it is asymptotic, just as the WKB expansion is. 19 Thus, 
for fixed N and 7, inclusion of additional coefficients in 
the expansion will eventually worsen the result. We can 
see this in Table II, where the error of our semiclassical 
result, |T SC — T|, at 7 = 1 is smaller than the error in 

the quadratic coefficients, iTsc 2 -* — T^\, and thus cannot 
be explained in terms of its approximation to that (or 
any higher) coefficient. On the other hand, the asymp- 

TABLE II: Coefficients of 7-expansion in Eq. (49) of the exact 
and semiclassical kinetic energy, T (i) and T$ , and the values 
T and T BC at 7 = 1 for v(x) — —10 sin 2 (ttx), where < x < 1. 



N 


T (0) 


T m 


T (2) 


rp(2) 
J SC 


T 




1 


2.31 


2.05 


0.614 


0.900 


5.07 


5.02 


2 


13.4 


9.78 


1.69 


1.62 


24.9 


24.7 



totic expansion with just the first few terms in Eq. (49) 
becomes accurate very rapidly as TV increases. Compare 
the relative error of the quadratic coefficients of about 
50% for N = 1 w ith roughly 4% for N = 2. The error 
dropped by about an order of magnitude as the number 
of particles increases by one. To understand why that is 
so, we next consider the A^-depcndence of each contribu- 
tion, where N is now the number of particles at 7 = 1. 
As N — > 00, the box must appear flat. Evaluating the 
local approximation on the flat box density yields: 



T, 



loc r 



TT 2 N 3 

6L 2 



1 



9 



(flat) 



(52) 

fixing the first three coeffcients with the values above, 
and the rest to vanish. The corrections to this flat limit 
can only involve powers of l/N, which we can either de- 
rive or find numerically. 
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Since the leading term is given by TF theory, if we 
expand the potential in the box in a power series around 
its average value, we find 



7T 2 iV 3 



6i 2 



3Sv 2 
(7m) 4 



where n = N/L and 



5v 2 = / dx(v(x) -vf/L 
Jo 



(53) 



(54) 



with v the ave rage of the potential over the well. For 
our shape, 5v 2 = D/8, yielding a TF value of 2.31, as 
in the figure. More importantly, we see that the leading 
correction to the flat box result is 0(1/N A ). Similarly, 
we find by fitting, that 



3ir 2 N 2 
16L 2 



1 



a 

N 2 



b_ 

iV 4 



(55) 



and is given exactly by the semiclassical approximation. 
For the specific choice of potential v(x) = — 10 sin 2 (7r:r) 
the coefficients are a = 0.38 and b = —0.26. Finally, 



T (2) 



16L 2 



1 



c 

iV 4 



(56) 



but the coefficient c is not given correctly by the semiclas- 
sical approximation. For v(x) = — 10 sin 2 (7ric) the exact 
value is c = 0.42, whereas the semiclassical approxima- 
tion gives about half that value. Thus, all corrections to 
the flat results vanish rapidly as N increases, and the er- 
rors of the first few semiclassical terms in the expansion 
become much smaller, leading to a much more accurate 
value at 7 = 1. 



TABLE III: Local approximation to the kinetic energy eval- 
uated on the TF, semiclassical, and exact density, and the 
kinetic energy from direct integration of the semiclassical and 
exact KED, all relative to the flat box value for N particles in 
a single-well potential v(x) — —10 sin 2 nx, < x < 1. The er- 
rors of our semiclassical result with respect to the exact result 
are denoted by Asc. 







AT loc 


[n\ 




AT 


N 


TF 


sc 


exact 


Asc 


sc 


exact 


Asc 


1 


-1.8 


0.96 


0.82 


0.14 


0.09 


0.13 


-0.04 


2 


-8.3 


0.89 


0.92 


-0.03 


0.08 


0.24 


-0.16 


4 


-32 


0.77 


0.78 


-0.01 


0.09 


0.14 


-0.05 


6 


-70 


0.72 


0.73 


-0.01 


0.07 


0.09 


-0.02 


8 


-123 


0.70 


0.70 


0.00 


0.06 


0.07 


-0.01 



In Table III, we list the various kinetic energies as func- 
tions of N for our well. Because the errors vanish so 
rapidly, we subtract the energies of the uniform system, 
as in: 



AT = T -T 



unif 



(57) 



and likewise for AT [n] . These differences could also be 
thought of as the change in energy due to turning on the 
well in the bottom of the box, analogous to the change 
in energy when atoms form a molecule. We see that our 
approximations become very accurate very quickly, and 
converge as 1/N 2 . 

The quantum correction yields a density that is not 
normalized. This is because the requirement in Eq. (23) 
that the phase vanishes at both x — and at x = L 
is used to determine ^ sc , not simple normalization. Of 
course, the error vanishes rapidly as N — > 00; for N = 1, 
it is AN = 4 x 10~ 2 and for N = 2, AiV = 6 x 10~ 4 . 
One can easily imagine schemes that patch this failure 
up, but we prefer to leave it as a measure of the overall 
error in the approximation. 

Since our formulas reduce to the exact results for a 
uniform potential, more generally, they should preserve 
these good features for a slowly-varying potential. We 
have applied our density formula to many examples, and 
almost always found it to be remarkably accurate. This 
is because of its excellent formal properties, and because 
we capture the leading correction to the LPA in a well- 
defined (albeit asymptotic) series. Most importantly, it 
appears that the conditions of application, fi sc above v 
everywhere, imply that these leading corrections always 
improve over the dominant contribution. 



B. Uniform convergence 

While the most important aspect of our work is the 
recovery of the leading asymptotic corrections to TF for 
the energies, the detailed spatial dependence is also im- 
portant for understanding how this is achieved, and also 
for understanding the strengths and weaknesses of this 
approach. 

Our semiclassical approximations are exact in the case 
of a uniform potential, where they yield the simple for- 
mulas: 



, unif 



(x) = N 1 



.(27ri\V)\ 



2A^ sin nx 



) 



(58) 



and Eq. (45) with / = l/sin(7T2;), w = sin(27rA^x), and 
N = (iV/7 + 1/2). These offer some insight into the 
nature of the expansion. 

Consider first the density. For any finite value of x, 
the oscillating contribution shrinks and oscillates more 
rapidly as N — > 00. Thus, we can expand the smooth 
part, the prefactor of the oscillating contribution, and 
the phase of the oscillation, in powers of 1/N, which is 
linear in 7 for small 7. On the other hand, for Nx fixed, 
one can again expand the density for large N: 



n 7 (y) = N 



1 - 



sin27ry\ ny sin(27ry) 



2ny 



12N 2 



(59) 



where y = Nx. The first term is precisely the profile of 
a semi-infinite box at the surface. This series is very ill- 
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behaved for large y, except for the lowest-order term. 
Similar comments apply to the kinetic energy density 
only more so, as several contributions diverge for small 
Nx. Thus there are two distinct regions and limits within 
the well, the interior and the edges. 

In what follows we illustrate that the error of our semi- 
classical approximations converges uniformly as 7 — > 0. 
We define: 



0.02 




-0.02 



An sc (x) = n sc (x) - n(x) 



(60) 



as the error in the semiclassical density, and likewise for 
the kinetic energy density. We pick a generic single- well 

0.002 1 , , , , 1 



-0.04 
-0.06 
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FIG. 7: Same as Fig. 5, for kinetic energy density. 
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FIG. 5: Fractional error in density for v(x) = —5 sin 2 (ti) 
only shown in interior (9 f /tt > 0.7 in left half). 
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FIG. 6: Fractional error in density close to the edge for v(x) = 
-5 sin 2 (the), 6>f/tt < 0.7. 

potential v(x) = — 5sin 2 (7ra;) sufficiently close to flat so 
that we are in a regime dominated by the asymptotic be- 
havior. For illustrative purposes we increase the extend 
of the edge-region by choosing [3 — 0.7. The fractional 



FIG. 8: Same as Fig. 6, for kinetic energy density. 



error of the density in the interior is shown in Fig. 5. 
The error converges uniformly throughout the interior as 
7 — > 0, being 0(7). As shown in Fig. 6, the fractional er- 
ror close to the edge of the box also converges uniformly, 
being 0(7) but noticeably larger. The convergence for 
the KED is shown in Figs. 7 and 8, and has the same 
features as the density, but is much larger. 



C. Phase oscillations 

We also check that the quantum oscillations of our 
semiclassical formulas can be extracted from the exact 
results as 7 — > 0. For a fixed point x we look at the 
difference between the exact result and the smooth term, 
multiplied by the prcfactor appearing in our formula for 
the quantum oscillations: 

<i 7 (x) = 2T F k F (x) sin a(x) An s , 7 (a;). (61) 
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If our results are correct, as 7 — > 0, this becomes a simple 
function of S F (x)/j, for any values ofxand7, specifically 
— sin(2S F (x)/j). The same analysis applies to the kinetic 




1 2 3 4 5 

2S f (x)/(ytc) 



above v max . This is the worst qualitative breakdown of 
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FIG. 9: Leading correction to the kinetic energy density for 
7 = 1, 0.5, 0.1 and v(x) = — 10 sin 2 (71-2;). 



FIG. 10: Exact and approximate ground-state densities for 
v(x) = —27 sin 2 (71-2;), where < x < 1. The lowest eigenvalue 
is eo = —16.3 and fj& c — 0.08. The position of the turning 
points is indicated by dashed lines. 



energy density, where we define: 



41^1^) 



(62) 



In Fig. 9, we plot g 7 (x) for several values of 7, as 
a function of 29 f (x)/tt, finding results converging to 
— sm.26 w (x), as predicted by the leading term of Eq. (42). 



D. Evanescent regions 

The only condition on the applicability of our approx- 
imations is that /j sc > v(x) for all x, But fi sc is between 
the highest occupied and lowest unoccupied level, and so 
for many well-depths, it can be the case that the HOMO 
has turning points while the condition is still valid. The 
starkest example is for N = 1, since beyond those turning 
points, the only occupied orbital is evanescent. Yet our 
approximations can still be applied, even though they 
contain only trigonometric functions of the phase, and 
no decaying exponentials, and still yield highly accurate 
results. In Figs. 1 and 3 we show results for a well depth 
of 12, for which the lowest eigenvalue is eo = —4.27 and 
/i sc = 5.52, and the turning points are located at around 
x = 0.2 and x = 0.8. Eventually the quadratic approach 
of the semiclassical density near the wall of the box mim- 
ics the exponential decay of the true density. Even the 
KED truncated by our method, only misses the nega- 
tive contribution, which largely cancels the error in the 
interior. The results for this well remain remarkably ac- 
curate. As the semiclassical chemical potential fi sc ap- 
proaches u max the validity of our approximation breaks 
down. We simulate such a situation in Figs. 10 and 11 by 
choosing a well depth of 27, such that fj, sc is only slightly 




FIG. 11: Exact and approximate ground-state KEDs for 
v(x) — —27 sin 2 (irx), where < x < 1. The lowest eigen- 
value is eo = —16.3 and /j, sc — 0.08. The position of the 
turning points is indicated by dashed lines. 

our approximation, yielding the errors shown in Tab. IV. 
But even here, errors are < 20%. 



CONSEQUENCES FOR DENSITY 
FUNCTIONAL THEORY 



This work has been confined to one-dimensional non- 
interacting particles confined by hard walls. In this sec- 
tion, we discuss in detail the ramifications for density 
functional theory in the real world of atoms, molecules, 
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TABLE IV: Exact and approximate quantities for one particle 
in a single well v(x) = — Dsin 2 (irx), where < x < 1. 



Energy levels 


£l 


£2 








D = 12 


-4.27 


13.6 


0.04 


5.52 




D = 27 


-16.3 


5.47 


-6.75 


0.08 




Kinetic energy 


T 


T ioc [n] 




exact 


sc 


TF 


exact 


sc 


D = 12 


5.13 


5.18 


2.66 


5.12 


5.33 


D = 27 


5.74 


7.63 


4.80 


6.42 


8.47 



and solids. We divide the discussion in two: Thomas- 
Fermi theory and Kohn-Sham theory. 

We begin with Thomas-Fermi theory and its exten- 
sions. This was the original density functional theory 
(DFT) and continues to be used in many fields of physics. 
TF theory became obsolescent for electronic structure 
calculations with Kohn-Sham work, but there has been 
a recent resurgence of interest in orbital- free DFT, with 
the hope of treating systems of much greater size than is 
presently possible with Kohn-Sham calculations. To do 
this, all that is needed is an accurate approximation to 
the non-interacting kinetic energy as a functional of the 
density. The original approximation using uniform gas 
inputs, is simply the 3D analog of our ID local approxi- 
mation used here. Thus if our methods could be gener- 
alized to apply to the general 3D case, it would produce 
an orbital-free theory. 

Perhaps the most important result of this study is to 
highlight an alternative path. Instead of trying to find 
density functional approximations, we have derived the 
leading corrections in terms of the potential, a perfectly 
valid alternative variable to the density. 37 If general for- 
mulas (or algorithms) could be found for finding accurate 
approximations to n[i; s ](r) and i s [u s ](r), where the sub- 
script s denotes non-interacting, one could use them to 
avoid solving the Kohn-Sham equations and evaluating 
any orbitals. At each step in the iteration, one finds v s (r), 
the Kohn-Sham potential, using some standard XC func- 
tional, and uses this to generate a new density. When 
self-consistency is reached, the kinetic energy is evalu- 
ated and the many-body energy is found in the usual 
way. A recent study 38 shows that this procedure is cor- 
rect once both approximations are derived from the same 
approximate Green's function, as ours have been. 

Even before such generalizations have been found, we 
have been able to use results here to deduce information 
on the 3D kinetic energy functional. In Id, our results 
show that the leading corrections to the asymptotic ex- 
pansion of the kinetic energy in powers of 1 /N are not de- 
termined by the gradient expansion for any finite system, 
but instead are given by the quantum corrections pro- 
ducing quantum oscillations. Ref. 39 is a careful study of 
the asymptotic expansion for the 3D kinetic energy, and 
showed how generalizing the gradient expansion to ensure 
recovery of the asymptotic expansion greatly improved 
total energies over the gradient expansion, but worsened 



other energies, such as those of jcllium surfaces. This 
reflects the difficulty in attempting to capture different 
physical limits with simple density-functional approxima- 
tions. Even our simple results cannot be easily encoded 
in a density functional approximation, but are both sim- 
ple and (relatively) physically transparent as potential 
functionals. 

Almost all modern electronic structure calculations 
are performed within the Kohn-Sham formalism, which 
provides a set of self-consistent non-interacting single- 
particle equations which reproduce the exact single- 
particle density of the interacting system. In these, the 
non-interacting kinetic energy is treated exactly, and only 
a small contribution to the total energy, the XC energy, 
is approximated as a functional of the density. This con- 
tribution is determined by the Coulomb repulsion, and 
so is a many-body effect. 

So, do we learn anything from studying our little toy 
problems? The answer is definitely yes. Our toy is 
perhaps the simplest possible system in which one can 
meaningfully approximate a Schrodingcr equation with 
its density functional analog, and make a local density 
approximation. So we learn in what limits this becomes 
relatively exact, and how to find the leading corrections. 
We learn the nature of these corrections (asymptotic) and 
how there are mutliple length scales in the system. While 
the details of these lessons depend on the functional we 
are approximating, some general features of functionals 
and their approximations can be guessed at, and highly 
useful analogies can be made. 

For example, there are many ways to understand 
why Kohn-Sham calculations are far more accurate than 
Thomas-Fermi type calculations, and our analysis pro- 
duces one more. A KS calculation, by virtue of its or- 
bitals, produces an incredibly accurate density, and we 
have seen how local-type approximations are much more 
powerful on accurate densities than on self-consistent 
densities. Thus not only does a Kohn-Sham calculation 
approximate only a small fraction of the total energy (the 
XC piece), but even that part is much more accurately 
given by a local approximation by virture of the accurate 
density. 

The insight based on the semiclassical analysis of func- 
tionals has already led to significant development in the 
Kohn-Sham XC functional. Schwinger demonstrated 40 ' 41 
that LDA exchange becomes relatively exact for large Z 
neutrals. Analysis of the large- Z behavior of modern 
exchange GGAs 42,43 shows that the most popular func- 
tionals all recover (to within about 20%) the leading cor- 
rections to the LDA asymptotic behavior of exchange for 
atoms. On the other hand, the gradient expansion ap- 
proximation, based on the slowly-varying limit, does not, 
being too small by almost exactly a factor of 2. This 
is entirely analogous to our problem, in which the local 
approximation recovers the exact dominant term, and 
a decent approximation to the next correction, but the 
gradient expansion worsens that agreement. Since such 
functionals are tested on exchange energy of atoms, and 
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these cannot be accurate without accurate asymptotic 
values, this is vital for recovering accurate thermochem- 
istry, which requires accurate atomic energies. On the 
other hand, bond lengths depend only on small variations 
in the energy when the bond distances is varied slightly 
around its equilibrium value, and so do not require accu- 
rate energies of isolated atoms. These can be improved 
upon over regular GGAs by restoring the true gradient 
expansion and ignoring the asymptotic limit. The recent 
PBEsol functional 28 does exactly this for exchange and 
produces better lattice parameters for many solids. 

VI. SUMMARY 

We have presented a fuller and more precise account 
of the results originally shown in Ref. 18. Kohn and 
Sham 14 produced asymptotic expansions for the interior, 
exterior, and turning point regions of the density. In 
Ref. 18, we presented a uniform approximation for the 
interior region, but only an asymptotic approximation 
for the kinetic energy density. Here, by analyzing the 
breakdown of the method for the boundary regions, we 
have produced a (nearly) uniform approximation to the 
kinetic energy. 

This work was supported by NSF under grant number 
CHE-0809859. We also acknowledge support from the 
KITP under grant number PHY05-51164. 



Thus, 

n sc (x) = f — g sc {x, e) = n s (x) + n osc (x) , (A5) 

where C(/i sc ) is a contour enclosing all poles of occupied 
states determined by /z sc . 

First we evaluate the density coming from the smooth 
term. In the limit L — > oo this is dominated by exp[— iO], 
simplifying the integral to 

n.(x) = ~<f tt~~~t i (A6) 

which, evaluated on the real axis, yields: 




Appendix A: Derivation of the semiclassical density 
and kinetic energy density in the complex plane 

The semiclassical corrections were derived from a con- 
tour integral over the semiclassical Green's function in 
Ref. 18, and we give a fuller account here. The method 
is well-described in Ref. 14 but we go beyond the aims 
there, since we require our solution to be uniformly 
asymptotic, not just producing the correct asymptotic 
expansion in the interior and we extract also the kinetic 
energy density. We are also treating box boundary con- 
ditions, rather than the turning points discussed there. 
Begin with the diagonal Green's function 



2cf> L (x) (/>r{x) 
W(e) 



(Al) 



where W(e) = <f>L{x)d4>n(x) / dx — 4>n(x)d(j)L(x)/dx is 
the Wronskian, and approximate the two independent 
solutions 4>l(x) and 4>r(x) via the WKB wavefunctions 
satisfying the boundary conditions: 



^ KB (x) 



wx[e(x)]/y/k(xj, 



^ KB (l) = Wx[0(L-x)]/y/k(x), 



(A2) 
(A3) 



yielding 

. cos0 — cos [26 (x) — Q] . . , , 
g sc {x,e) = . = g B {x,e)+g OBC {x,e) 



k(x) sin 



(A4) 



FIG. 12: Contour of integration in the complex e-plane. 

Then, we evaluate the oscillating term of Eq. (A4). 
We pick a contour C(/i sc ) as shown in Fig. 12, i.e., a 
vertical line along e = /i sc + z£ connected to a semi- 
circle, which encloses all poles of N lower-lying energy- 
eigenvalues ejy, . . . , ei. In the classical continuum limit 
fi sc >> C, allowing us to expand all quantities in the in- 
tegrand in powers of the imaginary part of the energy, 





1 



k((i sc + iQ, x) 
6(fj, sc + i(, x) 



1 - 



,(x) + iC,r F (x) + 



(A8) 
(A9) 



Keeping terms up to first order in £, employing the semi- 
classical quantization condition for the given boundary 
conditions in Eq. (23) with j = N + 1/2, and substitut- 
ing u = 4T F £, we obtain the result in Eq. (35). Note that 
the additional term of 1/2 in the quantization condition 
relative to Ref. 14 is due to the Mazlov index for a hard 
wall being 0, rather than 1/4 at a real turning point. 

Next, we provide some details of the lengthy derivation 
of the kinetic energy density of non-interacting, same- 
spin fermions: 



^sc (^0 



-A [e - v(x)]g sc (x, e) = t s (x) + t osc (x) . 

2lTl 

(A10) 
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In analogy to the derivation of the semiclassical density 
we first evaluate the smooth term yielding 



fcf(z) 

()7T 



(All) 



The subdominant piece of smooth term gives another 
contribution: 



de k(x, e) exp i0(e) 
27ri 4 sin 6(e) 



(A12) 



is evaluated on the contour C(/i sc ) as shown in Fig. 12. 
Hence, all quantities in the integrand are expanded in 
powers of £. In particular, we define s = — 2iT(Q, where 
r (C) = Jo dx [ y / 2[u, sc (x) + <] - y/2n sc (x)], express the 
C-expansion of all quantities in terms of s(C), and trun- 
cate its expansion after the quadratic term. Note that 
we approximate terms like [Jq dx/k F (x)]/T F by l/k%(x). 
This amounts to the same as neglecting the derivatives 
of Kj Then we obtain 



ds 



2itk F {x)T F l J z — e s 
where z = exp [2i0 F ] , which yields 

s { ' 24fc F (x)T F 2 ' 



(A13) 



(A14) 



where the total contribution of the smooth term is the 
sum of ti an d ti 2 \ agreement with Eq. (39). Then, 
we evaluate the oscillating term, which is integrated also 
along the contour C(/i sc ) in Fig. 12. 



We write the cosine of the oscillating piece as a 
weighted sum of exponential functions and demonstrate 
the derivation for the term that has the positive sign in 
the exponential function. We call this term to S \. 

As before we expand all quantities in powers of £. 
Here, we define q = — 41T(£), express the ^-expansion 
of all quantities in the integrand by q(C), and truncate 
its expansion after the quadratic term. Then we integrate 
by aid of the polygamma functions of order n, 44 defined 
as V (n) (0 = (-1)" +1 J™dq q n exp {-lq)/[l - exp (-?)], 
yielding 



k F (x) 



8ttT f 



y + l 



1 



16ttA: f (x)T f 2 
1 

' 1287rT|fc|(x) 



HI) 

*"'(f)-*'"(^ 



y + l 



(2) 

Similarly, the other term to SC integrates to the result in 
Eq. (A15) with y — > —y + l, where y = ot F (x)/n. The 
particular combination of the polygamma functions in 

£osc and i2c yields 



sin20 F (x) 
cos20 F (x) 
sin26»(Ai^) 



^osc (*^) 



fc F (x)sin20 F (x) 7r cosa F (a;) cos20 F (a;) 



4T F sina F (x) 4fc F (x)T 2 sin 2 a F (x 
7r 2 sin26» F (a;) /l 1 
8T^kf(x)sma F (x) \2 sin 2 a F (x) 



(A16) 



Finally, the sum of the smooth and oscillating pieces yield 
the result of Eq. (45). 
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